{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pricing European options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a two-step binomial tree. A non-dividend paying stock price starts at \\\\$50, and in each of the two time steps, the stock may go up by 20 percent or go down by 20 percent. We suppose that the risk-free rate is 5 percent per annum and the time to maturity T is 2 years. We would like to find the value of an European put option with a strike price K of $52."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "r = 0.05\n",
    "T = 2\n",
    "t = T/2\n",
    "u = 1.2\n",
    "d = 0.8\n",
    "\n",
    "q = (math.exp(r*t)-d)/(u-d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q is 0.6281777409400603\n"
     ]
    }
   ],
   "source": [
    "print('q is', q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing the StockOption base class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "\"\"\" \n",
    "Stores common attributes of a stock option \n",
    "\"\"\"\n",
    "class StockOption(object):\n",
    "    def __init__(\n",
    "        self, S0, K, r=0.05, T=1, N=2, pu=0, pd=0, \n",
    "        div=0, sigma=0, is_put=False, is_am=False):\n",
    "        \"\"\"\n",
    "        Initialize the stock option base class.\n",
    "        Defaults to European call unless specified.\n",
    "\n",
    "        :param S0: initial stock price\n",
    "        :param K: strike price\n",
    "        :param r: risk-free interest rate\n",
    "        :param T: time to maturity\n",
    "        :param N: number of time steps\n",
    "        :param pu: probability at up state\n",
    "        :param pd: probability at down state\n",
    "        :param div: Dividend yield\n",
    "        :param is_put: True for a put option,\n",
    "                False for a call option\n",
    "        :param is_am: True for an American option,\n",
    "                False for a European option\n",
    "        \"\"\"\n",
    "        self.S0 = S0\n",
    "        self.K = K\n",
    "        self.r = r\n",
    "        self.T = T\n",
    "        self.N = max(1, N)\n",
    "        self.STs = [] # Declare the stock prices tree\n",
    "\n",
    "        \"\"\" Optional parameters used by derived classes \"\"\"\n",
    "        self.pu, self.pd = pu, pd\n",
    "        self.div = div\n",
    "        self.sigma = sigma\n",
    "        self.is_call = not is_put\n",
    "        self.is_european = not is_am\n",
    "\n",
    "    @property\n",
    "    def dt(self):\n",
    "        \"\"\" Single time step, in years \"\"\"\n",
    "        return self.T/float(self.N)\n",
    "\n",
    "    @property\n",
    "    def df(self):\n",
    "        \"\"\" The discount factor \"\"\"\n",
    "        return math.exp(-(self.r-self.div)*self.dt)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A class for European options using a binomial tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "    import math\n",
    "    import numpy as np\n",
    "    from decimal import Decimal\n",
    "\n",
    "    \"\"\" \n",
    "    Price a European option by the binomial tree model \n",
    "    \"\"\"\n",
    "    class BinomialEuropeanOption(StockOption):\n",
    "\n",
    "        def setup_parameters(self):\n",
    "            # Required calculations for the model\n",
    "            self.M = self.N+1  # Number of terminal nodes of tree\n",
    "            self.u = 1+self.pu  # Expected value in the up state\n",
    "            self.d = 1-self.pd  # Expected value in the down state\n",
    "            self.qu = (math.exp(\n",
    "                (self.r-self.div)*self.dt)-self.d)/(self.u-self.d)\n",
    "            self.qd = 1-self.qu\n",
    "\n",
    "        def init_stock_price_tree(self):\n",
    "            # Initialize terminal price nodes to zeros\n",
    "            self.STs = np.zeros(self.M)\n",
    "\n",
    "            # Calculate expected stock prices for each node\n",
    "            for i in range(self.M):\n",
    "                self.STs[i] = self.S0 * \\\n",
    "                    (self.u**(self.N-i)) * (self.d**i)\n",
    "\n",
    "        def init_payoffs_tree(self):\n",
    "            \"\"\"\n",
    "            Returns the payoffs when the option \n",
    "            expires at terminal nodes\n",
    "            \"\"\" \n",
    "            if self.is_call:\n",
    "                return np.maximum(0, self.STs-self.K)\n",
    "            else:\n",
    "                return np.maximum(0, self.K-self.STs)\n",
    "\n",
    "        def traverse_tree(self, payoffs):\n",
    "            \"\"\"\n",
    "            Starting from the time the option expires, traverse\n",
    "            backwards and calculate discounted payoffs at each node\n",
    "            \"\"\"\n",
    "            for i in range(self.N):\n",
    "                payoffs = (payoffs[:-1]*self.qu + \n",
    "                           payoffs[1:]*self.qd)*self.df\n",
    "\n",
    "            return payoffs\n",
    "\n",
    "        def begin_tree_traversal(self):\n",
    "            payoffs = self.init_payoffs_tree()\n",
    "            return self.traverse_tree(payoffs)\n",
    "\n",
    "        def price(self):\n",
    "            \"\"\" Entry point of the pricing implementation \"\"\"\n",
    "            self.setup_parameters()\n",
    "            self.init_stock_price_tree()\n",
    "            payoffs = self.begin_tree_traversal()\n",
    "            \n",
    "            # Option value converges to first node\n",
    "            return payoffs[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "eu_option = BinomialEuropeanOption(\n",
    "    50, 52, r=0.05, T=2, N=2, pu=0.2, pd=0.2, is_put=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "European put option price is: 4.1926542806038585\n"
     ]
    }
   ],
   "source": [
    "print('European put option price is:', eu_option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A class for American options using a binomial tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "\n",
    "\"\"\" \n",
    "Price a European or American option by the binomial tree \n",
    "\"\"\"\n",
    "class BinomialTreeOption(StockOption):\n",
    "\n",
    "    def setup_parameters(self):\n",
    "        self.u = 1+self.pu  # Expected value in the up state\n",
    "        self.d = 1-self.pd  # Expected value in the down state\n",
    "        self.qu = (math.exp(\n",
    "            (self.r-self.div)*self.dt)-self.d)/(self.u-self.d)\n",
    "        self.qd = 1-self.qu\n",
    "\n",
    "    def init_stock_price_tree(self):\n",
    "        # Initialize a 2D tree at T=0\n",
    "        self.STs = [np.array([self.S0])]\n",
    "\n",
    "        # Simulate the possible stock prices path\n",
    "        for i in range(self.N):\n",
    "            prev_branches = self.STs[-1]\n",
    "            st = np.concatenate(\n",
    "                (prev_branches*self.u, \n",
    "                 [prev_branches[-1]*self.d]))\n",
    "            self.STs.append(st) # Add nodes at each time step\n",
    "\n",
    "    def init_payoffs_tree(self):\n",
    "        if self.is_call:\n",
    "            return np.maximum(0, self.STs[self.N]-self.K)\n",
    "        else:\n",
    "            return np.maximum(0, self.K-self.STs[self.N])\n",
    "\n",
    "    def check_early_exercise(self, payoffs, node):\n",
    "        if self.is_call:\n",
    "            return np.maximum(payoffs, self.STs[node] - self.K)\n",
    "        else:\n",
    "            return np.maximum(payoffs, self.K - self.STs[node])\n",
    "\n",
    "    def traverse_tree(self, payoffs):\n",
    "        for i in reversed(range(self.N)):\n",
    "            # The payoffs from NOT exercising the option\n",
    "            payoffs = (payoffs[:-1]*self.qu + \n",
    "                       payoffs[1:]*self.qd)*self.df\n",
    "\n",
    "            # Payoffs from exercising, for American options\n",
    "            if not self.is_european:\n",
    "                payoffs = self.check_early_exercise(payoffs,i)\n",
    "\n",
    "        return payoffs\n",
    "\n",
    "    def begin_tree_traversal(self):\n",
    "        payoffs = self.init_payoffs_tree()\n",
    "        return self.traverse_tree(payoffs)\n",
    "\n",
    "    def price(self):\n",
    "        \"\"\"  The pricing implementation \"\"\"\n",
    "        self.setup_parameters()\n",
    "        self.init_stock_price_tree()\n",
    "        payoffs = self.begin_tree_traversal()\n",
    "        return payoffs[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "am_option = BinomialTreeOption(50, 52, \n",
    "    r=0.05, T=2, N=2, pu=0.2, pd=0.2, is_put=True, is_am=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "American put option price is: 5.089632474198373\n"
     ]
    }
   ],
   "source": [
    "print(\"American put option price is:\", am_option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Cox-Ross-Rubinstein model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A class for Cox-Ross-Rubinstein binomial tree option pricing model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "\"\"\" \n",
    "Price an option by the binomial CRR model \n",
    "\"\"\"\n",
    "class BinomialCRROption(BinomialTreeOption):\n",
    "    def setup_parameters(self):\n",
    "        self.u = math.exp(self.sigma * math.sqrt(self.dt))\n",
    "        self.d = 1./self.u\n",
    "        self.qu = (math.exp((self.r-self.div)*self.dt) - \n",
    "                   self.d)/(self.u-self.d)\n",
    "        self.qd = 1-self.qu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider again the two-step binomial tree. The non-dividend paying stock has a current price of \\\\$50 and a volatility of 30 percent. Suppose that the risk-free rate is 5 percent per annum and the time to maturity T is 2 years. We would like to find the value of an European put option with a strike price K of \\\\$52 by the CRR model:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "eu_option = BinomialCRROption(\n",
    "    50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "European put: 6.245708445206436\n"
     ]
    }
   ],
   "source": [
    "print('European put:', eu_option.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "am_option = BinomialCRROption(\n",
    "    50, 52, r=0.05, T=2, N=2, \n",
    "    sigma=0.3, is_put=True, is_am=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "American put option price is: 7.428401902704834\n"
     ]
    }
   ],
   "source": [
    "print('American put option price is:', am_option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using a Leisen-Reimer (LR) tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A class for LR binomial tree option pricing model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "\"\"\" \n",
    "Price an option by the Leisen-Reimer tree\n",
    "\"\"\"\n",
    "class BinomialLROption(BinomialTreeOption):\n",
    "\n",
    "    def setup_parameters(self):\n",
    "        odd_N = self.N if (self.N%2 == 0) else (self.N+1)\n",
    "        d1 = (math.log(self.S0/self.K) +\n",
    "              ((self.r-self.div) +\n",
    "               (self.sigma**2)/2.)*self.T)/\\\n",
    "            (self.sigma*math.sqrt(self.T))\n",
    "        d2 = (math.log(self.S0/self.K) +\n",
    "              ((self.r-self.div) -\n",
    "               (self.sigma**2)/2.)*self.T)/\\\n",
    "            (self.sigma * math.sqrt(self.T))\n",
    "\n",
    "        pbar = self.pp_2_inversion(d1, odd_N)\n",
    "        self.p = self.pp_2_inversion(d2, odd_N)\n",
    "        self.u = 1/self.df * pbar/self.p\n",
    "        self.d = (1/self.df-self.p*self.u)/(1-self.p)\n",
    "        self.qu = self.p\n",
    "        self.qd = 1-self.p\n",
    "\n",
    "    def pp_2_inversion(self, z, n):\n",
    "        return .5 + math.copysign(1, z)*\\\n",
    "            math.sqrt(.25 - .25*\n",
    "                math.exp(\n",
    "                    -((z/(n+1./3.+.1/(n+1)))**2.)*(n+1./6.)\n",
    "                )\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "eu_option = BinomialLROption(\n",
    "    50, 52, r=0.05, T=2, N=4, sigma=0.3, is_put=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "European put: 5.878650106601964\n"
     ]
    }
   ],
   "source": [
    "print(\"European put:\", eu_option.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "am_option = BinomialLROption(50, 52, \n",
    "    r=0.05, T=2, N=4, sigma=0.3, is_put=True, is_am=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "American put: 6.763641952939979\n"
     ]
    }
   ],
   "source": [
    "print('American put:', am_option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Greeks for free"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A class for greeks with the LR binomial tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\"\"\" \n",
    "Compute option price, delta and gamma by the LR tree \n",
    "\"\"\"\n",
    "class BinomialLRWithGreeks(BinomialLROption):\n",
    "\n",
    "    def new_stock_price_tree(self):\n",
    "        \"\"\"\n",
    "        Creates an additional layer of nodes to our\n",
    "        original stock price tree\n",
    "        \"\"\"\n",
    "        self.STs = [np.array([self.S0*self.u/self.d,\n",
    "                              self.S0,\n",
    "                              self.S0*self.d/self.u])]\n",
    "\n",
    "        for i in range(self.N):\n",
    "            prev_branches = self.STs[-1]\n",
    "            st = np.concatenate((prev_branches*self.u,\n",
    "                                 [prev_branches[-1]*self.d]))\n",
    "            self.STs.append(st)\n",
    "\n",
    "    def price(self):\n",
    "        self.setup_parameters()\n",
    "        self.new_stock_price_tree()\n",
    "        payoffs = self.begin_tree_traversal()\n",
    "\n",
    "        # Option value is now in the middle node at t=0\n",
    "        option_value = payoffs[len(payoffs)//2]\n",
    "\n",
    "        payoff_up = payoffs[0]\n",
    "        payoff_down = payoffs[-1]\n",
    "        S_up = self.STs[0][0]\n",
    "        S_down = self.STs[0][-1]\n",
    "        dS_up = S_up - self.S0\n",
    "        dS_down = self.S0 - S_down\n",
    "\n",
    "        # Calculate delta value\n",
    "        dS = S_up - S_down\n",
    "        dV = payoff_up - payoff_down\n",
    "        delta = dV/dS\n",
    "\n",
    "        # calculate gamma value\n",
    "        gamma = ((payoff_up-option_value)/dS_up - \n",
    "                 (option_value-payoff_down)/dS_down) / \\\n",
    "            ((self.S0+S_up)/2. - (self.S0+S_down)/2.)\n",
    "\n",
    "        return option_value, delta, gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "eu_call = BinomialLRWithGreeks(50, 52, r=0.05, T=2, N=300, sigma=0.3)\n",
    "results = eu_call.price()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "European call values\n",
      "Price: 9.69546807138366\n",
      "Delta: 0.6392477816643529\n",
      "Gamma: 0.01764795890533088\n"
     ]
    }
   ],
   "source": [
    "print('European call values')\n",
    "print('Price: %s\\nDelta: %s\\nGamma: %s' % results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "eu_put = BinomialLRWithGreeks(\n",
    "    50, 52, r=0.05, T=2, N=300, sigma=0.3, is_put=True)\n",
    "results = eu_put.price()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "European put values\n",
      "Price: 6.747013809252746\n",
      "Delta: -0.3607522183356649\n",
      "Gamma: 0.0176479589053312\n"
     ]
    }
   ],
   "source": [
    "print('European put values')\n",
    "print('Price: %s\\nDelta: %s\\nGamma: %s' % results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trinomial trees in options pricing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  A class for the trinomial tree options pricing model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "\n",
    "class TrinomialTreeOption(BinomialTreeOption):\n",
    "\n",
    "    def setup_parameters(self):\n",
    "        \"\"\" Required calculations for the model \"\"\"\n",
    "        self.u = math.exp(self.sigma*math.sqrt(2.*self.dt))\n",
    "        self.d = 1/self.u\n",
    "        self.m = 1\n",
    "        self.qu = ((math.exp((self.r-self.div) *\n",
    "                             self.dt/2.) -\n",
    "                    math.exp(-self.sigma *\n",
    "                             math.sqrt(self.dt/2.))) /\n",
    "                   (math.exp(self.sigma *\n",
    "                             math.sqrt(self.dt/2.)) -\n",
    "                    math.exp(-self.sigma *\n",
    "                             math.sqrt(self.dt/2.))))**2\n",
    "        self.qd = ((math.exp(self.sigma *\n",
    "                             math.sqrt(self.dt/2.)) -\n",
    "                    math.exp((self.r-self.div) *\n",
    "                             self.dt/2.)) /\n",
    "                   (math.exp(self.sigma *\n",
    "                             math.sqrt(self.dt/2.)) -\n",
    "                    math.exp(-self.sigma *\n",
    "                             math.sqrt(self.dt/2.))))**2.\n",
    "\n",
    "        self.qm = 1 - self.qu - self.qd\n",
    "\n",
    "    def init_stock_price_tree(self):\n",
    "        # Initialize a 2D tree at t=0\n",
    "        self.STs = [np.array([self.S0])]\n",
    "\n",
    "        for i in range(self.N):\n",
    "            prev_nodes = self.STs[-1]\n",
    "            self.ST = np.concatenate(\n",
    "                (prev_nodes*self.u, [prev_nodes[-1]*self.m,\n",
    "                                     prev_nodes[-1]*self.d]))\n",
    "            self.STs.append(self.ST)\n",
    "\n",
    "    def traverse_tree(self, payoffs):\n",
    "        # Traverse the tree backwards \n",
    "        for i in reversed(range(self.N)):\n",
    "            payoffs = (payoffs[:-2] * self.qu +\n",
    "                       payoffs[1:-1] * self.qm +\n",
    "                       payoffs[2:] * self.qd) * self.df\n",
    "\n",
    "            if not self.is_european:\n",
    "                payoffs = self.check_early_exercise(payoffs,i)\n",
    "\n",
    "        return payoffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "eu_put = TrinomialTreeOption(\n",
    "    50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "European put: 6.573565269142496\n"
     ]
    }
   ],
   "source": [
    "print('European put:', eu_put.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "am_option = TrinomialTreeOption(50, 52, \n",
    "    r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "American put: 7.161349217272585\n"
     ]
    }
   ],
   "source": [
    "print('American put:', am_option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lattices in options pricing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A class for CRR binomial lattice option pricing model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "class BinomialCRRLattice(BinomialCRROption):\n",
    "\n",
    "    def setup_parameters(self):\n",
    "        super(BinomialCRRLattice, self).setup_parameters()\n",
    "        self.M = 2*self.N + 1\n",
    "\n",
    "    def init_stock_price_tree(self):\n",
    "        self.STs = np.zeros(self.M)\n",
    "        self.STs[0] = self.S0 * self.u**self.N\n",
    "\n",
    "        for i in range(self.M)[1:]:\n",
    "            self.STs[i] = self.STs[i-1]*self.d\n",
    "\n",
    "    def init_payoffs_tree(self):\n",
    "        odd_nodes = self.STs[::2]  # Take odd nodes only\n",
    "        if self.is_call:\n",
    "            return np.maximum(0, odd_nodes-self.K)\n",
    "        else:\n",
    "            return np.maximum(0, self.K-odd_nodes)\n",
    "\n",
    "    def check_early_exercise(self, payoffs, node):\n",
    "        self.STs = self.STs[1:-1]  # Shorten ends of the list\n",
    "        odd_STs = self.STs[::2]  # Take odd nodes only\n",
    "        if self.is_call:\n",
    "            return np.maximum(payoffs, odd_STs-self.K)\n",
    "        else:\n",
    "            return np.maximum(payoffs, self.K-odd_STs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "eu_option = BinomialCRRLattice(\n",
    "    50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "European put: 6.245708445206432\n"
     ]
    }
   ],
   "source": [
    "print('European put:', eu_option.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "am_option = BinomialCRRLattice(50, 52, \n",
    "    r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "American put: 7.428401902704828\n"
     ]
    }
   ],
   "source": [
    "print(\"American put:\", am_option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing the TrinomialLattice class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\"\"\" \n",
    "Price an option by the trinomial lattice \n",
    "\"\"\"\n",
    "class TrinomialLattice(TrinomialTreeOption):\n",
    "\n",
    "    def setup_parameters(self):\n",
    "        super(TrinomialLattice, self).setup_parameters()\n",
    "        self.M = 2*self.N + 1\n",
    "\n",
    "    def init_stock_price_tree(self):\n",
    "        self.STs = np.zeros(self.M)\n",
    "        self.STs[0] = self.S0 * self.u**self.N\n",
    "\n",
    "        for i in range(self.M)[1:]:\n",
    "            self.STs[i] = self.STs[i-1]*self.d\n",
    "\n",
    "    def init_payoffs_tree(self):\n",
    "        if self.is_call:\n",
    "            return np.maximum(0, self.STs-self.K)\n",
    "        else:\n",
    "            return np.maximum(0, self.K-self.STs)\n",
    "\n",
    "    def check_early_exercise(self, payoffs, node):\n",
    "        self.STs = self.STs[1:-1]  # Shorten ends of the list\n",
    "        if self.is_call:\n",
    "            return np.maximum(payoffs, self.STs-self.K)\n",
    "        else:\n",
    "            return np.maximum(payoffs, self.K-self.STs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "European put: 6.573565269142496\n"
     ]
    }
   ],
   "source": [
    "eu_option = TrinomialLattice(\n",
    "    50, 52, r=0.05, T=2, N=2, sigma=0.3, is_put=True)\n",
    "print('European put:', eu_option.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "American put: 7.161349217272585\n"
     ]
    }
   ],
   "source": [
    "am_option = TrinomialLattice(50, 52, \n",
    "    r=0.05, T=2, N=2, sigma=0.3, is_put=True, is_am=True)\n",
    "print('American put:', am_option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite differences in options pricing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing the Finite Differences base class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "from abc import ABC, abstractmethod\n",
    "import numpy as np\n",
    "\n",
    "\"\"\" \n",
    "Base class for sharing \n",
    "attributes and functions of FD \n",
    "\"\"\"\n",
    "class FiniteDifferences(object):\n",
    "\n",
    "    def __init__(\n",
    "        self, S0, K, r=0.05, T=1, \n",
    "        sigma=0, Smax=1, M=1, N=1, is_put=False\n",
    "    ):\n",
    "        self.S0 = S0\n",
    "        self.K = K\n",
    "        self.r = r\n",
    "        self.T = T\n",
    "        self.sigma = sigma\n",
    "        self.Smax = Smax\n",
    "        self.M, self.N = M, N\n",
    "        self.is_call = not is_put\n",
    "\n",
    "        self.i_values = np.arange(self.M)\n",
    "        self.j_values = np.arange(self.N)\n",
    "        self.grid = np.zeros(shape=(self.M+1, self.N+1))\n",
    "        self.boundary_conds = np.linspace(0, Smax, self.M+1)\n",
    "\n",
    "    @property\n",
    "    def dS(self):\n",
    "        return self.Smax/float(self.M)\n",
    "\n",
    "    @property\n",
    "    def dt(self):\n",
    "        return self.T/float(self.N)\n",
    "\n",
    "    @abstractmethod\n",
    "    def setup_boundary_conditions(self):\n",
    "        raise NotImplementedError('Implementation required!')\n",
    "\n",
    "    @abstractmethod\n",
    "    def setup_coefficients(self):\n",
    "        raise NotImplementedError('Implementation required!')\n",
    "\n",
    "    @abstractmethod\n",
    "    def traverse_grid(self):\n",
    "        \"\"\"  \n",
    "        Iterate the grid backwards in time \n",
    "        \"\"\"\n",
    "        raise NotImplementedError('Implementation required!')\n",
    "\n",
    "    @abstractmethod\n",
    "    def interpolate(self):\n",
    "        \"\"\"\n",
    "        Use piecewise linear interpolation on the initial\n",
    "        grid column to get the closest price at S0.\n",
    "        \"\"\"\n",
    "        return np.interp(\n",
    "            self.S0, self.boundary_conds, self.grid[:,0])\n",
    "\n",
    "    def price(self):\n",
    "        self.setup_boundary_conditions()\n",
    "        self.setup_coefficients()\n",
    "        self.traverse_grid()\n",
    "        return self.interpolate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A class for pricing European options using explicit method of finite differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\"\"\" \n",
    "Explicit method of Finite Differences \n",
    "\"\"\"\n",
    "class FDExplicitEu(FiniteDifferences):\n",
    "\n",
    "    def setup_boundary_conditions(self):\n",
    "        if self.is_call:\n",
    "            self.grid[:,-1] = np.maximum(\n",
    "                0, self.boundary_conds - self.K)\n",
    "            self.grid[-1,:-1] = (self.Smax-self.K) * \\\n",
    "                np.exp(-self.r*self.dt*(self.N-self.j_values))\n",
    "        else:\n",
    "            self.grid[:,-1] = np.maximum(\n",
    "                0, self.K-self.boundary_conds)\n",
    "            self.grid[0,:-1] = (self.K-self.Smax) * \\\n",
    "                np.exp(-self.r*self.dt*(self.N-self.j_values))\n",
    "\n",
    "    def setup_coefficients(self):\n",
    "        self.a = 0.5*self.dt*((self.sigma**2) *\n",
    "                              (self.i_values**2) -\n",
    "                              self.r*self.i_values)\n",
    "        self.b = 1 - self.dt*((self.sigma**2) *\n",
    "                              (self.i_values**2) +\n",
    "                              self.r)\n",
    "        self.c = 0.5*self.dt*((self.sigma**2) *\n",
    "                              (self.i_values**2) +\n",
    "                              self.r*self.i_values)\n",
    "\n",
    "    def traverse_grid(self):\n",
    "        for j in reversed(self.j_values):\n",
    "            for i in range(self.M)[2:]:\n",
    "                self.grid[i,j] = \\\n",
    "                    self.a[i]*self.grid[i-1,j+1] +\\\n",
    "                    self.b[i]*self.grid[i,j+1] + \\\n",
    "                    self.c[i]*self.grid[i+1,j+1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the example of an European put option. The underlying stock price is $50 with a volatility of 40 percent. The strike price of the put option is $50 with an expiration time of 5 months. The risk-free rate is 10 percent.\n",
    "\n",
    "We can price the option using the explicit method with a Smax value of 100, an M value of 100, and a N value of 1,000:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.072882278148043\n"
     ]
    }
   ],
   "source": [
    "option = FDExplicitEu(50, 50, r=0.1, T=5./12., \n",
    "    sigma=0.4, Smax=100, M=100, N=1000, is_put=True)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens when other values of M and N are chosen improperly?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-8.109445694129245e+35\n"
     ]
    }
   ],
   "source": [
    "option = FDExplicitEu(50, 50, r=0.1, T=5./12., \n",
    "    sigma=0.4, Smax=100, M=80, N=100, is_put=True)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A class for pricing European options using the implicit method of finite differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as linalg\n",
    "\n",
    "\"\"\" \n",
    "Explicit method of Finite Differences \n",
    "\"\"\"\n",
    "class FDImplicitEu(FDExplicitEu):\n",
    "\n",
    "    def setup_coefficients(self):\n",
    "        self.a = 0.5*(self.r*self.dt*self.i_values -\n",
    "                      (self.sigma**2)*self.dt*\\\n",
    "                          (self.i_values**2))\n",
    "        self.b = 1 + \\\n",
    "                 (self.sigma**2)*self.dt*\\\n",
    "                    (self.i_values**2) + self.r*self.dt\n",
    "        self.c = -0.5*(self.r*self.dt*self.i_values +\n",
    "                       (self.sigma**2)*self.dt*\\\n",
    "                           (self.i_values**2))\n",
    "        self.coeffs = np.diag(self.a[2:self.M],-1) + \\\n",
    "                      np.diag(self.b[1:self.M]) + \\\n",
    "                      np.diag(self.c[1:self.M-1],1)\n",
    "\n",
    "    def traverse_grid(self):\n",
    "        \"\"\" Solve using linear systems of equations \"\"\"\n",
    "        P, L, U = linalg.lu(self.coeffs)\n",
    "        aux = np.zeros(self.M-1)\n",
    "\n",
    "        for j in reversed(range(self.N)):\n",
    "            aux[0] = np.dot(-self.a[1], self.grid[0, j])\n",
    "            x1 = linalg.solve(L, self.grid[1:self.M, j+1]+aux)\n",
    "            x2 = linalg.solve(U, x1)\n",
    "            self.grid[1:self.M, j] = x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.071594188049893\n"
     ]
    }
   ],
   "source": [
    "option = FDImplicitEu(50, 50, r=0.1, T=5./12., \n",
    "    sigma=0.4, Smax=100, M=100, N=1000, is_put=True)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.063684691731647\n"
     ]
    }
   ],
   "source": [
    "option = FDImplicitEu(50, 50, r=0.1, T=5./12., \n",
    "    sigma=0.4, Smax=100, M=80, N=100, is_put=True)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A class for pricing European options using the Crank-Nicolson method of finite differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as linalg\n",
    "\n",
    "\"\"\" \n",
    "Crank-Nicolson method of Finite Differences \n",
    "\"\"\"\n",
    "class FDCnEu(FDExplicitEu):\n",
    "\n",
    "    def setup_coefficients(self):\n",
    "        self.alpha = 0.25*self.dt*(\n",
    "            (self.sigma**2)*(self.i_values**2) - \\\n",
    "            self.r*self.i_values)\n",
    "        self.beta = -self.dt*0.5*(\n",
    "            (self.sigma**2)*(self.i_values**2) + self.r)\n",
    "        self.gamma = 0.25*self.dt*(\n",
    "            (self.sigma**2)*(self.i_values**2) +\n",
    "            self.r*self.i_values)\n",
    "        self.M1 = -np.diag(self.alpha[2:self.M], -1) + \\\n",
    "                  np.diag(1-self.beta[1:self.M]) - \\\n",
    "                  np.diag(self.gamma[1:self.M-1], 1)\n",
    "        self.M2 = np.diag(self.alpha[2:self.M], -1) + \\\n",
    "                  np.diag(1+self.beta[1:self.M]) + \\\n",
    "                  np.diag(self.gamma[1:self.M-1], 1)\n",
    "\n",
    "    def traverse_grid(self):\n",
    "        \"\"\" Solve using linear systems of equations \"\"\"\n",
    "        P, L, U = linalg.lu(self.M1)\n",
    "\n",
    "        for j in reversed(range(self.N)):\n",
    "            x1 = linalg.solve(\n",
    "                L, np.dot(self.M2, self.grid[1:self.M, j+1]))\n",
    "            x2 = linalg.solve(U, x1)\n",
    "            self.grid[1:self.M, j] = x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.072238354486825\n"
     ]
    }
   ],
   "source": [
    "option = FDCnEu(50, 50, r=0.1, T=5./12.,\n",
    "    sigma=0.4, Smax=100, M=100, N=1000, is_put=True)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.070145703042843\n"
     ]
    }
   ],
   "source": [
    "option = FDCnEu(50, 50, r=0.1, T=5./12., \n",
    "    sigma=0.4, Smax=100, M=80, N=100, is_put=True)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pricing exotic barrier options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A down-and-out option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\"\"\"\n",
    "Price a down-and-out option by the Crank-Nicolson\n",
    "method of finite differences.\n",
    "\"\"\"\n",
    "class FDCnDo(FDCnEu):\n",
    "\n",
    "    def __init__(\n",
    "        self, S0, K, r=0.05, T=1, sigma=0, \n",
    "        Sbarrier=0, Smax=1, M=1, N=1, is_put=False\n",
    "    ):\n",
    "        super(FDCnDo, self).__init__(\n",
    "            S0, K, r=r, T=T, sigma=sigma,\n",
    "            Smax=Smax, M=M, N=N, is_put=is_put\n",
    "        )\n",
    "        self.barrier = Sbarrier\n",
    "        self.boundary_conds = \\\n",
    "            np.linspace(Sbarrier, Smax, M+1)\n",
    "        self.i_values = self.boundary_conds/self.dS\n",
    "\n",
    "    @property\n",
    "    def dS(self):\n",
    "        return (self.Smax-self.barrier)/float(self.M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.491560552934787\n"
     ]
    }
   ],
   "source": [
    "option = FDCnDo(50, 50, r=0.1, T=5./12., \n",
    "    sigma=0.4, Sbarrier=40, Smax=100, M=120, N=500)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5413635028954452\n"
     ]
    }
   ],
   "source": [
    "option = FDCnDo(50, 50, r=0.1, T=5./12., sigma=0.4, \n",
    "    Sbarrier=40, Smax=100, M=120, N=500, is_put=True)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pricing American options with finite differences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A class for pricing American options using the Crank-Nicolson method of finite differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "\"\"\" \n",
    "Price an American option by the Crank-Nicolson method \n",
    "\"\"\"\n",
    "class FDCnAm(FDCnEu):\n",
    "\n",
    "    def __init__(self, S0, K, r=0.05, T=1, sigma=0, \n",
    "            Smax=1, M=1, N=1, omega=1, tol=0, is_put=False):\n",
    "        super(FDCnAm, self).__init__(S0, K, r=r, T=T, \n",
    "            sigma=sigma, Smax=Smax, M=M, N=N, is_put=is_put)\n",
    "        self.omega = omega\n",
    "        self.tol = tol\n",
    "        self.i_values = np.arange(self.M+1)\n",
    "        self.j_values = np.arange(self.N+1)\n",
    "\n",
    "    def setup_boundary_conditions(self):\n",
    "        if self.is_call:\n",
    "            self.payoffs = np.maximum(0, \n",
    "                self.boundary_conds[1:self.M]-self.K)\n",
    "        else:\n",
    "            self.payoffs = np.maximum(0, \n",
    "                self.K-self.boundary_conds[1:self.M])\n",
    "\n",
    "        self.past_values = self.payoffs\n",
    "        self.boundary_values = self.K * np.exp(\n",
    "                -self.r*self.dt*(self.N-self.j_values))\n",
    "        \n",
    "    def calculate_payoff_start_boundary(self, rhs, old_values):\n",
    "        payoff = old_values[0] + \\\n",
    "            self.omega/(1-self.beta[1]) * \\\n",
    "                (rhs[0] - \\\n",
    "                 (1-self.beta[1])*old_values[0] + \\\n",
    "                 self.gamma[1]*old_values[1])\n",
    "\n",
    "        return max(self.payoffs[0], payoff)      \n",
    "    \n",
    "    def calculate_payoff_end_boundary(self, rhs, old_values, new_values):\n",
    "        payoff = old_values[-1] + \\\n",
    "            self.omega/(1-self.beta[-2]) * \\\n",
    "                (rhs[-1] + \\\n",
    "                 self.alpha[-2]*new_values[-2] - \\\n",
    "                 (1-self.beta[-2])*old_values[-1])\n",
    "\n",
    "        return max(self.payoffs[-1], payoff)\n",
    "    \n",
    "    def calculate_payoff(self, k, rhs, old_values, new_values):\n",
    "        payoff = old_values[k] + \\\n",
    "            self.omega/(1-self.beta[k+1]) * \\\n",
    "                (rhs[k] + \\\n",
    "                 self.alpha[k+1]*new_values[k-1] - \\\n",
    "                 (1-self.beta[k+1])*old_values[k] + \\\n",
    "                 self.gamma[k+1]*old_values[k+1])\n",
    "\n",
    "        return max(self.payoffs[k], payoff)\n",
    "\n",
    "    def traverse_grid(self):\n",
    "        \"\"\" Solve using linear systems of equations \"\"\"\n",
    "        aux = np.zeros(self.M-1)\n",
    "        new_values = np.zeros(self.M-1)\n",
    "\n",
    "        for j in reversed(range(self.N)):\n",
    "            aux[0] = self.alpha[1]*(self.boundary_values[j] +\n",
    "                                    self.boundary_values[j+1])\n",
    "            rhs = np.dot(self.M2, self.past_values) + aux\n",
    "            old_values = np.copy(self.past_values)\n",
    "            error = sys.float_info.max\n",
    "\n",
    "            while self.tol < error:\n",
    "                new_values[0] = \\\n",
    "                    self.calculate_payoff_start_boundary(\n",
    "                        rhs, old_values)            \n",
    "\n",
    "                for k in range(self.M-2)[1:]:\n",
    "                    new_values[k] = \\\n",
    "                        self.calculate_payoff(\n",
    "                            k, rhs, old_values, new_values)                    \n",
    "\n",
    "                new_values[-1] = \\\n",
    "                    self.calculate_payoff_end_boundary(\n",
    "                        rhs, old_values, new_values)\n",
    "\n",
    "                error = np.linalg.norm(new_values-old_values)\n",
    "                old_values = np.copy(new_values)\n",
    "\n",
    "                self.past_values = np.copy(new_values)\n",
    "\n",
    "        self.values = np.concatenate(\n",
    "            ([self.boundary_values[0]], new_values, [0]))\n",
    "\n",
    "    def interpolate(self):\n",
    "        # Use linear interpolation on final values as 1D array\n",
    "        return np.interp(\n",
    "            self.S0, self.boundary_conds, self.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's price an American call-and-put option with an underlying asset price of 50 and volatility of 40 percent, a strike price of 50, a risk-free rate of 10 percent, and\n",
    "an expiration date of 5 months. We choose a <i>S<sub>max</sub></i> value of 100, <i>M</i> as 100, <i>N</i> as 42, an\n",
    "`omega` parameter value of 1.2, and a tolerance value of 0.001:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.108682815392217\n"
     ]
    }
   ],
   "source": [
    "option = FDCnAm(50, 50, r=0.1, T=5./12., \n",
    "    sigma=0.4, Smax=100, M=100, N=42, omega=1.2, tol=0.001)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.277764229383736\n"
     ]
    }
   ],
   "source": [
    "option = FDCnAm(50, 50, r=0.1, T=5./12., sigma=0.4, Smax=100, \n",
    "    M=100, N=42, omega=1.2, tol=0.001, is_put=True)\n",
    "print(option.price())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implied volatilities of AAPL American put option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Get implied volatilities from a Leisen-Reimer binomial\n",
    "tree using the bisection method as the numerical procedure.\n",
    "\"\"\"\n",
    "class ImpliedVolatilityModel(object):\n",
    "\n",
    "    def __init__(self, S0, r=0.05, T=1, div=0, \n",
    "                 N=1, is_put=False):\n",
    "        self.S0 = S0\n",
    "        self.r = r\n",
    "        self.T = T\n",
    "        self.div = div\n",
    "        self.N = N\n",
    "        self.is_put = is_put\n",
    "\n",
    "    def option_valuation(self, K, sigma):\n",
    "        \"\"\" Use the binomial Leisen-Reimer tree \"\"\"\n",
    "        lr_option = BinomialLROption(\n",
    "            self.S0, K, r=self.r, T=self.T, N=self.N, \n",
    "            sigma=sigma, div=self.div, is_put=self.is_put\n",
    "        )\n",
    "        return lr_option.price()\n",
    "\n",
    "    def get_implied_volatilities(self, Ks, opt_prices):\n",
    "        impvols = []\n",
    "        for i in range(len(strikes)):\n",
    "            # Bind f(sigma) for use by the bisection method\n",
    "            f = lambda sigma: \\\n",
    "                self.option_valuation(Ks[i], sigma)-\\\n",
    "                opt_prices[i]\n",
    "            impv = bisection(f, 0.01, 0.99, 0.0001, 100)[0]\n",
    "            impvols.append(impv)\n",
    "\n",
    "        return impvols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bisection(f, a, b, tol=0.1, maxiter=10):\n",
    "    \"\"\"\n",
    "    :param f: The function to solve\n",
    "    :param a: The x-axis value where f(a)<0\n",
    "    :param b: The x-axis value where f(b)>0\n",
    "    :param tol: The precision of the solution\n",
    "    :param maxiter: Maximum number of iterations\n",
    "    :return: The x-axis value of the root,\n",
    "                number of iterations used\n",
    "    \"\"\"\n",
    "    c = (a+b)*0.5  # Declare c as the midpoint ab\n",
    "    n = 1  # Start with 1 iteration\n",
    "    while n <= maxiter:\n",
    "        c = (a+b)*0.5\n",
    "        if f(c) == 0 or abs(a-b)*0.5 < tol:\n",
    "            # Root is found or is very close\n",
    "            return c, n\n",
    "\n",
    "        n += 1\n",
    "        if f(c) < 0:\n",
    "            a = c\n",
    "        else:\n",
    "            b = c\n",
    "                \n",
    "    return c, n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "strikes = [75, 80, 85, 90, 92.5, 95, 97.5, \n",
    "           100, 105, 110, 115, 120, 125]\n",
    "put_prices = [0.16, 0.32, 0.6, 1.22, 1.77, 2.54, 3.55, \n",
    "              4.8, 7.75, 11.8, 15.96, 20.75, 25.81]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = ImpliedVolatilityModel(\n",
    "    99.62, r=0.0248, T=78/365., div=0.0182, N=77, is_put=True)\n",
    "impvols_put = model.get_implied_volatilities(strikes, put_prices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd4VHX2+PH3SaOEGkJvCV1AqSKCiGBDUbGsHfuu61pXXcuuruv61S02FNuqP+uuveuKC6ggEhDpXUIKEIqQAoFQUs/vj3ujY5hkJpCZO5mc1/PMk5lbz52ZzLmfcu9HVBVjjDGmJjFeB2CMMSbyWbIwxhgTkCULY4wxAVmyMMYYE5AlC2OMMQFZsjDGGBOQJQsTFiJyv4j8x33eTUSKRCT2ELaTIiIqInEhiPFKEZl7GOt/ISJX+NuWe7w9alh3tYiccKj7Diff46xm/qUiMiNE+47I98n3+x2tLFnUkojMFpGdItKomvn3uz9mI6pMv1JEyt0fjd0iskxEznDnnSAim2ux/wPudvJE5EMR6RjkuioivWqYf1g/lsFS1U2q2kxVy+tyuyIyXUQe8DN9koj8WJcJxt+Pg6qepqqv+VvePd4sd91XReTBKvMHqOrsuoovlGo6Tnf+G6p6Soj2fUjvk88Jiu9DReR2n2VuEpFs9/9zkYgcV6fB13OWLGpBRFKAMYACZ/mZL8BlQAHg78xrvqo2A1oBLwHvikjSIYRyo7udPu62phzCNqLRq8Bl7ufg6zLgDVUtC39IDUsoSnx1wecEpZn7v3MkUAF8ACAixwD/AH4FtMT5//zoUEq/0cqSRe1cDnyH86PkLxmMAToBtwAXiUiCv42oagXwMtAEqLZqIhBVLcD5sg+En0odv66c71tSEJE57uTl7lnVhYG2LyIbROQOEVkhIntF5CURae9WQ+wRkS9FpLW7bGX10LUislVEtvmetVXZ7i+qkkSkpbvtbSKyRUQerPwnFZFYEXnULUVlARNrCPljIAnnc6jcV2vgDOB1n329LiK5IrJRRO4VEb//ByLypIjkuGeai0VkjDt9AvAn4EL3vVzuTv/F+19lWyoivUTkWuBS4E533c983uuT3OcxInK3iGSKSL6I/HRSISKNReQ/7vRdIrJQRNpXs89OIvKBe6zZInKzz7xpIvKYz+t3RORl9/mVIpImIk+JSKGI/CAiJ/os+9Nx+iw7RUQKgPurllDdY79ORNaLUyp/pjKhu5/vY+7nmy0iN0oN1YxV3qf73ffmdff7uFpEhvtbz4/LgTmqusF9nQKsVtXF6tzW4nUgGWhXTRypIvKNu9+Z7rK+898TpzRbKCJzRGSAO/1oEdnue3wicp6ILHOfjxCnVLPbXe7xII8n5CxZ1M7lwBvu41Q//6RXAJ8B77ivz/C3EfeL8mugCFh/qMGISDJwHrA00LKqerz7dJB7dvVOjSv87DzgZJxSzJnAFzg/lMk435+bqyw/DugNnALcXfmPHcBrQBnQCxjirlv5o/sbnPdxCDAc58zPL1XdD7yL8zlVugD4QVWXu6+fwjlz7AGMdZe9qppNLgQG4ySgN4H3RKSxqv4P+BvwjvteDgriGCtjfAHn+/Owu+6Zfha7GTjbja8TsBN4xp13hRt/V6ANcB2wv+oG3AT4GbAc6AycCPxeRE51F7kapxQ2XkQuBY7GOcmpdAyQhfM5/wX4UKovBVcu2w54qJplznD3MQjnM6mM4zfAaTjv81D3uGvjLOBtnBL2p8DTQa53Oc73rtIXQKyIHOOeqFwNLAN+rGb9N4HFOO/P/3HwyeMXOP8H7YAlOJ85qroQyMf5n6o0Gfi3+/xJ4ElVbQH0xPk+RwRLFkESp/6yO/Cuqi4GMoFLfOY3Bc4H3lTVUuB9Dv4CjRSRXThfwIuBc1S18BDCmepuZzmwDbjtELYRrKdUdbuqbgG+BRao6lJVLQY+wvkR9/VXVd2rqiuBV3COs1puwj0N+L273g6carWL3EUuAJ5Q1Ry3JPX3APG+BpwvIk3c1z/9KLg/AhcCf1TVPe5Z5WM41VQHUdX/qGq+qpap6mNAI6BvgP3Xhd8C96jqZvd9vh/4lXuSUYqTJHqparl7JrzbzzaOBtqq6gOqWuK2l7yI+76q6o84ieY1nB+oy1V1j8/6O3De91L3xGId1ZfqtqrqU+77dFDicv1DVXep6iZgFk5yAOfzfdI91p04VUG1MVdVp7ntX//GSUY1ckuI7XH+RyvtwSmlzwWKcRLktern5nki0g3n/f2zqhar6hycxPwTVX3Z/Y5Vfn6DRKSlO/s1nASBm4BPxUk+4Hy+vUQkWVWLVPW7YN6EcLBkEbwrgBmqmue+fpNfJoNzcM6Op7mv3wBOE5G2Pst8p6qtVDVZVUeq6peHGMvN7nY6q+qlqpp7iNsJxnaf5/v9vG5WZfkcn+cbcc6Ma9IdiAe2udUqu4Dn+bn438nPNqulqnOBXGCSOL2Pjubnf8RkIKHKNjbinHkfRERuF5G1blXCLpwz+mR/y9ax7jj15ZXvx1qgHOcH7t/AdOBtcar7HhaR+Gq20alyG+52/uRuo9J/gVhgnfu++dpS5Yeyps8yp5rpvnzP0Pfx8/em6ucbzLZq2m7j6qqwfFwBfKCqRT7Tfo1TmhiA8x2ZDPxXRPwdcydgp6ru9Zn203fKrVr7h1uNuBvY4M6q/O78BzhTRJrhJMtvVXWbO+8anFL8D24Vo9/aCS9EZGNUpHHPUi/AKaZWfjkbAa1EZJBbxXEFzj/ApsrqWJwfwYuBqWEKdS/Q1Od1hzDt11dX4Af3eTdga4Dlc3DO5JKraYDe5m6zUrcgYngdp0TRFyfBVya4PJwzt+7AGp/tbam6Affs8y6c6pvVqlohIjtxPldwOjkcqkDr5gBXq2paNfP/CvxVnA4X03DO+l/ys41sVe1dw34ewklEqSJysaq+5TOvs4iIT8LohlPN48/hvBfbgC4+r7tWt2BdcP+Xz8c5ufM1CPhMVdPd1/8TkW3AKH5ZAgEn5tYikuiTMLrx8/twCTAJOAknUbTEqUoUAFXdIiLz3RguA56r3LCqrgcudqsRzwXeF5E2VRKTJ6xkEZyzcc7s+uMUnwcDR+BUy1wuIpV1wmf4zB8E/BP/DeF+idN46fuo2qsnkGXAuSLSVJwustdUmb+dw2hQD9Kf3f0PwGkLqLFtxD2jmgE8JiItxGnc7SkiY91F3gVuFpEu4jRW3x1EDK/j/KP+Bp96abeq4l3gIRFpLiLdcarw/PWPb45TUswF4kTkPqCFz/ztQIpU0zgeQKDP4V9ujN0BRKStiExyn48TkSPdKrXdOMnPXxfk74HdInKXiDRxz3YHisjR7naOx/l8LncfT7nf40rtcN73eBE5H+f7Po269y5wi4h0FpFWOAk6lM4BduFUhflaCEwUkR7iqGynW1V1A6q6EViEk7AT3Cpq37an5jgnQPk4J29/8xPH68CdOL2yPqqcKCKTRaSt2wlmlzu5TruYHypLFsG5AnjF7X73Y+UDpzHtUpx/umWqOqPK/KnAUSIyMIh9dMap1vF99KxlnFOAEpwfo9dwG9V83A+85lZLXFDLbQfrGyAD+Ap4VFWDuTjrcpyi/xqcM7D3gcprR17EqXZZjtNQ+GGgjbltEfOARA4+G74JpwSWhVM//SZOz7SqpuM0UqbjVDEc4JdVJO+5f/NFZEmgmKp4Cejvfg4f+5n/pBv3DBHZg9MD7xh3Xgec92c3TqngG/wkOzcxnolz4pKNU6r6f0BLEWmB82N1o6pucaugXgJe8TlBWYDTQJuHUwL5larm1/I4g/EizsnCCpyOGtNwknSofiCvAF730xbxOk5D+Wyc93Yq8FtV/QH/LsH5TApw2jder7KtjTgl1jU4n19VH+FWN1YpNUwAVotIEc734CJVPRD00YWQ+Gm/MabW3CqRbCC+muokU0+IyJXAr1U17BelichpwL9UtXu49x1uIpKJk5AOte0yrKxkYYzxjFtFdrqIxLnVYH/Bp1omWonIeThtHF97HUuwrIHbGOMlwWmwfwen6vVz4D5PIwoxEZmN0/55mds2US9YNZQxxpiArBrKGGNMQFFTDZWcnKwpKSleh2GMMfXK4sWL81S1baDloiZZpKSksGjRIq/DMMaYekVEarwrQiWrhjLGGBOQJQtjjDEBWbIwxhgTkCULY4wxAVmyMMYYE5AlC2OMMQFZsjDGGBOQJQtjjKnHPli8mXcWbgr5fixZGGNMPfbit1l8sizQgJSHz5KFMcbUU3lFxfzw4x5G9wr90PCWLIwxpp6an+kMXjiqZ5uQ78uShTHG1FNpGXk0bxTHkZ1bhnxfliyMMaaeSsvMY2TPNsTFhv6nPKR7EJEJIrJORDJE5G4/868TkZUiskxE5opI/yrzu4lIkYj8IZRxGmNMfZNTsI+cgv2MDkMVFIQwWYhILPAMcBrOEIIXV00GwJuqeqSqDgYeBh6vMn8K8EWoYjTGmPoqLSMPICyN2xDaksUIIENVs1S1BHgbmOS7gKru9nmZiDOAOQAicjaQBawOYYzGGFMvzc3Io13zRvRq1yws+wtlsugM5Pi83uxO+wURuUFEMnFKFje70xKBu3AGcjfGGOOjokKZn5nP6F7JiEhY9hnKZOHvCPSgCarPqGpPnORwrzv5r8AUVS2qcQci14rIIhFZlJube9gBG2NMfbBu+x7y95aEpctspVAOq7oZ6OrzugtQ02WGbwPPuc+PAX4lIg8DrYAKETmgqk/7rqCqLwAvAAwfPvygRGSMMdEo3O0VENpksRDoLSKpwBbgIuAS3wVEpLeqrndfTgTWA6jqGJ9l7geKqiYKY4xpqNIy8uiRnEinVk3Cts+QJQtVLRORG4HpQCzwsqquFpEHgEWq+ilwo4icBJQCO4ErQhWPMcZEg9LyChZkF3Du0IOagEMqlCULVHUaMK3KtPt8nt8SxDbur/vIjDGmflqWs4t9JeWM7hm+KiiwK7iNMaZeScvIQwSODWPjNliyMMaYemVeRj4DO7WkVdOEsO7XkoUxxtQTe4vLWJqzk1G9wluqAEsWxhhTb3y/oYDScg17ewVYsjDGmHpjXkYeCbExHJ2SFPZ9W7Iwxph6Ii0jn6HdW9EkITbs+7ZkYYwx9UB+UTFrtu32pAoKLFkAsHbbbg6UlnsdhjHGVGt+ljOE6ujeliw8kZVbxMSp3/LS3GyvQzHGmGqlZeTTvFEcR4VhCFV/Gnyy6NG2GScd0Z5nZ2WwY/cBr8Mxxhi/5mXmcUyPpLAMoepPg08WAH86/QhKyit4dMY6r0MxxpiD5BTsY2P+vrDeZbYqSxZASnIiV45K4b3Fm1m1pdDrcIwx5hfmZYb/luRVWbJw3Ti+N62bJvB//12Dqg2NYYyJHGkZ+bRt3ojeYRpC1R9LFq6WTeK59eQ+LMguYPrqH70OxxhjAFBV5mXmMapnm7ANoeqPJQsfFx/dlT7tm/G3aT9QXGZdaY0x3lu3fQ95RSWeVkGBJYtfiIuN4c9n9GdTwT5eTdvgdTjGGENahnt9hSWLyDKmd1vG92vHU19nkFdU7HU4xpgGbl5GHiltmtI5jEOo+mPJwo8/nX4EB0rLeXxmutehGGMasNLyCr7Lyve8VAGWLPzq1a4Zlx3bnbe/38Tabbu9DscY00Ct2LyLvSXlliwi2S0n9qZ543ge/Ny60hpjvJGWke8Modoj/IMdVWXJohqtmiZw60m9ScvI56u1O7wOxxjTAM3NyKN/xxa0TgzvEKr+WLKowaUju9OzbSIPTVtLSVmF1+EYYxqQfSVlLN20k+MioAoKLFnUKD42hnsn9ic7by///m6j1+EYYxqQhRt2UlqujLJkUT+c0Lctx/dpy5NfplOwt8TrcIwxDcS8jDziY4WjU1p7HQpgySIgEeHeiUewt6ScJ760rrTGmPCYm5HH0G6taZoQ53UogCWLoPRp35xLRnTjjQWbWL99j9fhGGOi3M69Jc4QqhFSBQWWLIJ268l9aJoQy4Ofr/U6FGNMlJuflY8qjO7lfZfZSpYsgpSUmMAtJ/bmm/RcZq2zrrTGmNCZm5FHs0ZxHNWlldeh/MSSRS1cfmwKKW2a8tDnayktt660xpjQmJeRxzGpScR7NISqP5ETST2QEBfDPRP7k7GjiDcXbPI6HGNMFNq8cx8b8vdFTJfZSpYsaumkI9oxqmcbpnyZTuG+Uq/DMcZEmXk/3ZI8ctorwJJFrYkIfz6jP7v3l/LkV+u9DscYE2XSMvNIbpZA3/bNvQ7lFyxZHIIjOrbgwqO78fr8DWTmFnkdjjEmSjhDqOYzqmeyp0Oo+hPSZCEiE0RknYhkiMjdfuZfJyIrRWSZiMwVkf7u9JNFZLE7b7GIjA9lnIfitpP70Dg+lr9ZV1pjTB1Zv6OI3D3FEVcFBSFMFiISCzwDnAb0By6uTAY+3lTVI1V1MPAw8Lg7PQ84U1WPBK4A/h2qOA9V2+aNuHF8L776YQffrs/1OhxjTBSYuz4P8H4IVX9CWbIYAWSoapaqlgBvA5N8F1BV35GFEgF1py9V1a3u9NVAYxFpFMJYD8lVo1PoltSUB/+7ljLrSmuMOUzzMvPo3qYpXVo39TqUg4QyWXQGcnxeb3an/YKI3CAimTgli5v9bOc8YKmqRtyA2I3iYvnjaf1Yt30P7yzKCbyCMcZUo6y8gu+yChjVM/JKFRDaZOGvdeagIedU9RlV7QncBdz7iw2IDAD+CfzW7w5ErhWRRSKyKDfXm6qgCQM7MCI1icdmpLP7gHWlNcYcmuWbCykqLovI9goIbbLYDHT1ed0F2FrNsuBUU51d+UJEugAfAZeraqa/FVT1BVUdrqrD27ZtWwch156IcN8Z/dm5r4Snv87wJAZjTP03L8Npr2iIJYuFQG8RSRWRBOAi4FPfBUSkt8/LicB6d3or4HPgj6qaFsIY68TAzi351dAuvJKWzYa8vV6HY4yph9IynSFUkyJgCFV/apUsRCRGRFoEs6yqlgE3AtOBtcC7qrpaRB4QkbPcxW4UkdUisgy4DafnE+56vYA/u91ql4lIu9rEGm53nNqX+NgY/v6FdaU1xtTO/pJylmzcFbFVUAABR9UQkTeB64ByYDHQUkQeV9VHAq2rqtOAaVWm3efz/JZq1nsQeDDQ9iNJuxaNuWFcLx6Zvo55mXkRW5Q0xkSehRsKKCmviMgus5WCKVn0d7u4no3zw98NuCykUdVT1xyXSudWTXjwv2sprzioLd8YY/xKy3SGUB2RmuR1KNUKJlnEi0g8TrL4RFVL8dOryUDj+FjuPq0fa7bt5v3F1pXWGBOceRn5DOkaOUOo+hNMsnge2IBz0dwcEekO7K5xjQbsjKM6Mqx7ax6Znk5RcZnX4RhjItyufSWs2loY0VVQEESyUNWpqtpZVU9Xx0ZgXBhiq5cq70qbV1TMs7OsK60xpmbzMyNvCFV/AiYLEWkvIi+JyBfu6/783GvJ+DG4ayvOHdKZ/zc3m5yCfV6HY4yJYGmZeSQmxDKoa+QMoepPMNVQr+J0f+3kvk4Hfh+qgKLFHRP6EiPwj//94HUoxpgINi8jnxERNoSqP8FEl6yq7wIV8NP1E+UhjSoKdGzZhOvG9uTzFdtYuKHA63CMMRFo6679ZOXtjfj2CgguWewVkTa4PaBEZCRQGNKoosRvj+9Jx5aNeeCzNVRYV1pjTBVpGZF7S/KqgkkWt+HcpqOniKQBrwM3hTSqKNEkIZY7J/Rl5ZZCPlq6xetwjDERZl5mPm0SI28IVX+C6Q21BBgLjMK5++sAVV0R6sCixaRBnRnUtRUPT/+BvdaV1hjjUlXmZuQxqlcyMTGRNYSqP9Umi8qhTEXkXOAsoC/QBzjTnWaCEBPj3JV2++5inv/G781zjTENUEblEKo9I7vLbKWaLhccC3wNnOlnngIfhiSiKDSse2vOHNSJ5+dkceGIbnRu1cTrkIwxHqtP7RVQQ7JQ1b+4Tx9Q1WzfeSKSGtKootBdE/oyY/WPPPy/H3jyoiFeh2OM8djcjHy6JjWha1LkDaHqTzAN3B/4mfZ+XQcS7bq0bsq1x/fgk2VbWbJpp9fhGGM8VFZewYKsfI6rJ6UKqLnNop+InIdzS/JzfR5XAo3DFmEUuW5sT9o1b8QDn61B1brSGtNQrdxSyJ7isno1lEFNJYu+wBlAK5x2i8rHUOA3oQ8t+iQ2iuOOU/uyLGcXny6vaYRZY0w0m5eZD8CoetK4DTW3WXwCfCIix6rq/DDGFNXOG9qF1+Zv4J9f/MAp/TvQJCHW65CMMWE2d30eR3RsQZtmjbwOJWg1VUPd6T69RESmVn2EKb6o43SlHcDWwgO8+G2W1+EYY8LsQGk5izftrDddZivV1HW2cjDpReEIpCEZkZrE6Ud24LnZmVwwvCsdWloTkDENxaINOykpi+whVP2pqRrqM/fva+ELp+H442lH8OWaHTwyfR2PXTDI63CMMWEyNyOPuJjIHkLVn2qThYh8Rg3Dp6rqWSGJqIHomtSUq0an8MK3WVw3tge968G9YYwxh29eZh5DurUisVHkDqHqT03RPhq2KBqo68b25I0Fm5jyZTrPXjrM63CMMSFWuK+UlVsKuXl8b69DqbWaqqG+CWcgDVHrxASuHp3C1K8zWL21kAGdWnodkjEmhOZnVQ6hWr/aK6Dm3lDvun9XisiKqo/whRjdrhnTgxaN45gyM93rUIwxIZaWkUfThFgGR/gQqv7UVA11i/v3jHAE0lC1bBLPtcf34NEZ6SzL2VUvv0TGmOCkZeYxIjWJhLjIHkLVn2ojVtVt7tPrVXWj7wO4PjzhNQxXjk4lKTGBx2as8zoUY0yIbCvcT1buXkbXo1t8+AomvZ3sZ9ppdR1IQ9asURzXje3Bt+vz+D7bxus2JhqlZTi3+KiP7RVQc5vF70RkJdC3SntFNmBtFnXsspEptG3eiMdmrLObDBoTheZl5JGUmEC/DvWzm3xNJYs3cW4c+Cm/vJHgMFWdHIbYGpQmCbHccEJPFmQX/HSTMWNMdFBV0jLzOLZnm3oxhKo/NbVZFKrqBlW92G2n2I9zkV4zEekWtggbkItGdKNjy8ZWujAmymTm7mX77uJ6214BQbRZiMiZIrIeyAa+ATYAX4Q4rgapcXwsN43vzZJNu5i9LtfrcIwxdaRyCNX6NNhRVcE0cD8IjATSVTUVOBFIC2lUDdj5w7vQNakJj8200oUx0SItI48urZvQrU39GELVn2CSRamq5gMxIhKjqrOAwSGOq8GKj43hlhP7sGrLbqav3u51OMaYw1ReoczPyq/XVVAQXLLYJSLNgDnAGyLyJFAWzMZFZIKIrBORDBG528/869wrxJeJyFwR6e8z74/ueutE5NRgDyganD24Ez2SE5kyM52KCitdGFOfrdxSyJ4DZYzuHf3JYhJO4/atwP+ATJxeUTUSkVjgGZxrMvoDF/smA9ebqnqkqg4GHgYed9ftD1wEDAAmAM+622sQ4mJj+P3JfVi3fQ//Xbkt8ArGmIhV2V5Rn4ZQ9SdgslDVvaparqplqvqaqk51q6UCGQFkqGqWqpYAb+MkHt9t7/Z5mcjPt0SfBLytqsWqmg1kuNtrMM44siN92zfniS/TKSuv8DocY8whmpeZR78OzUmuR0Oo+lPTRXl7RGS3z2OP798gtt0ZyPF5vdmdVnU/N4hIJk7J4ubarBvNYmKEW0/uTVbuXj5ettXrcIwxh+BAaTkLN+yst1dt+6rpOovmqtrC59Hc928Q2/Z35clBFfCq+oyq9gTuAu6tzboicq2ILBKRRbm50dfV9NQBHRjQqQVPfpVOqZUujKl3Fm+sHEK1fldBQXBtFojIIBG50X0cFeS2NwNdfV53AWo6RX4bOLs266rqC6o6XFWHt23bNsiw6g8R4fZT+pBTsJ/3Fm32OhxjTC2l/TSEagNIFiJyC/AG0M59vCEiNwWx7YVAbxFJFZEEnAbrT6ts23e4qInAevf5p8BFItJIRFKB3sD3Qewz6ozr244h3Vrx9NfrKS4r9zocY0wtpGXmM6hrK5rVsyFU/QmmZHENcIyq3qeq9+FcoPebQCupahlwIzAdWAu8q6qrReQBEakcv/tGEVktIsuA24Ar3HVXA+8Ca3B6YN2gqg3yl1JEuP3kvmwtPMDb3+cEXsEYExEK95eycvOuqGivgJoHP6okgO8PdTn+2xQOoqrTgGlVpt3n8/yWg1b6ed5DwEPB7Cfaje7VhmNSk3h6VgYXDO9Kk4QG04vYmHrru6x8KhRG1/Mus5WCKVm8AiwQkftF5H7gO+ClkEZlfsFpu+hL7p5i/vPdRq/DMcYEYV5GHk3iYxnSrbXXodSJYK6zeBy4GigAdgJXqeoToQ7M/NKI1CTG9E7muW8yKSoO6gJ6Y4yH5mbU3yFU/anpOovPReRSEUlU1cXuxXhPqurScAZofnb7KX0p2FvCa/M2eB2KMaYGPxYeIDN3b1R0ma1UU8p7ATgD2CAi74jI2W6vJuORwV1bcWK/djz/TSaF+0u9DscYU415mZW3+IiOxm2o+aK8T1T1YqAb8CFOT6VNIvKyiPgbl9uEwa0n92H3gTJempvtdSjGmGqkZeTTumk8/TsGc/1y/RBMm8V+VX1HVc8BTgGG4HRnNR4Y2Lklpw3swMtzs9m5t8TrcIwxVagqaRl5jOqZXG+HUPUnmIvy2ovITSKSBnwMzACGhTwyU61bT+7D3pIynp+T5XUoxpgqsvL28uPuA4yKovYKqLmB+zci8jWwBOgD3KmqPVT1LlVdFrYIzUH6tG/OWYM68dq8DeTuKfY6HGOMj3nuLcnr+2BHVdVUshgF/APoqqo3qaoNpRpBbjmxNyXlFTw3O9PrUIwxPuZm5NG5VRO61+MhVP2pqYH7KlWdoap2u9MI1KNtM84d0pn/LNjIj4UHvA7HGIM7hGpmPqN7tUEketorIMi7zprIdPOJvVFVnp61PvDCxpiQW721kN0HyqLmflC+LFnUY12TmnLB8K68szCHnIJ9XodjTIOXluEMInpslNwPyldNDdxJNT3CGaSp3o3jeyEiPPW1lS6M8VpaRh592zenXfPGXodS52oqWSwGFrkLoEyMAAAeGUlEQVR/c4F0nPEmct1pJgJ0bNmES4/pxgdLtpCdt9frcIxpsJwhVAuirstspZoauFNVtQfOeBRnqmqyqrbBuQXIh+EK0AT2uxN6Eh8rPPllutehGNNgLdm0k+KyiqjrMlspmDaLo91xKQBQ1S+AsaELydRWu+aNuWJUCp8s38r67Xu8DseYBiktI4/YGOGYHtFZSx9MssgTkXtFJEVEuovIPUB+qAMztfPb43vSND6WKVa6MMYTaRn5DOrSkuaN470OJSSCSRYXA22Bj9xHW3eaiSBJiQlcc1wq01b+yOqthV6HY0yDUrivlBVRNISqP8HcSLDAHf50jKoOVdXfq2pBGGIztXTNmB60aBzHlJnWM8qYcHopLZsKhQkDO3gdSsgEcyPBUSKyBljjvh4kIs+GPDJTay2bxHPt8T34cu12luXs8jocYxqE/KJiXvo2i9OP7MCATi29DidkgqmGmgKcittOoarLgeNDGZQ5dFeOTqV103gen2ltF8aEw7OzM9lfWs5tJ/f1OpSQCuoKblXNqTKpPASxmDrQrFEc143tyZz0XBZusNpCY0Jp6679/Pu7jZw3tAu92jXzOpyQCiZZ5IjIKEBFJEFE/gCsDXFc5jBcfmwKyc0a8ej0daiq1+EYE7We/HI9KPz+5D5ehxJywSSL64AbgM7AZmCw+9pEqCYJsdwwricLsguYl2m9nI0JhczcIt5bnMOlI7vRuVUTr8MJuWB6Q+Wp6qWq2l5V26nqZFW1X6AId/GIbnRs2ZjHZljpwphQeHxGOo3jY7lhXC+vQwmLuOpmiMidqvqwiDwFHPRro6o3hzQyc1gax8dy4/he3PPRKmavy2Vcv3Zeh2RM1Fi1pZDPV27j5vG9SG7WyOtwwqKmkkVlu0TlzQSrPkyEO39YV7omNeGxmVa6MKYuPTJ9Ha2axvPr43t4HUrYVFuyUNXP3L+vhS8cU5cS4mK4eXxv7nh/BdNXb4/qC4aMCZcFWfl8k57Ln07vR4sovbWHPzVVQ32Gn+qnSqp6VkgiMnXqnCGdeW52JlNmpnNK//bExETXUI/GhJOq8vD0dbRv0YjLj03xOpywqjZZAI+GLQoTMnGxMdxyUm9ueXsZn6/cxpmDOnkdkjH11tc/7GDxxp08dM5AGsfHeh1OWNU0nsU3lQ9gPrATKADmu9NMPXHmUZ3o074ZU75Mp6y8wutwjKmXKiqUR6avI6WNM5xxQxPMvaEmApnAVOBpIENETgt1YKbuxMQIt53ch6zcvXyybKvX4RhTL322Yis//LiHW0/uQ3xsUDe/iCrBHPFjwDhVPUFVxwLjcO4XZeqRUwd0YECnFjz51XpKrXRhTK2Ullfw+Mx0jujYgjOPaphVucEkix2qmuHzOgvYEczGRWSCiKwTkQwRudvP/NtEZI2IrBCRr0Sku8+8h0VktYisFZGpImIts4dBRLj9lD5sKtjH+4s3ex2OMfXKOwtz2Ji/jztO7dNgO4kEkyxWi8g0EblSRK4APgMWisi5InJudSuJSCzwDHAa0B+4WET6V1lsKTBcVY8C3gcedtcdBYwGjgIGAkdjQ7ketnF92zG4ayue+mo9xWV2L0hjgnGgtJypX61nePfWjOvbcC9uDSZZNAa24/xYnwDkAknAmcAZNaw3AshQ1SxVLQHeBib5LqCqs1R1n/vyO6BL5Sx3vwlAIyDejcEcBhHhD6f0ZWvhAd7+vuqNhI0x/rw2bwM79hRz54R+NOQKjpq6zgKgqlcd4rY7A76/SJuBY2pY/hrgC3ef80VkFrANEOBpVbU73daB0b3aMCI1iadnZXDB8K40SWhY3f+MqY3dB0p57ptMxvZpy4jUJK/D8VQwvaFSReRxEflQRD6tfASxbX8p2O9FfiIyGRgOPOK+7gUcgVPS6AyMF5GDBlwSkWtFZJGILMrNzQ0iJFNZusjdU8zfv7D8a0xNXpyTxa59pdxxanQPbBSMgCUL4GPgJZy2itp0o9kM+HZG7gIc1G9TRE4C7gHGqmqxO/kc4DtVLXKX+QIYCczxXVdVXwBeABg+fLjd/ChII1KT+PVxqfy/udkcnZJkF+oZ40funmJempvNxKM6MrBz9A6XGqxg2iwOqOpUt33B90K9QBYCvd2SSQJwEfCLEomIDAGeB85SVd8eVpuAsSISJyLxOO0ldhpch+46rR9Du7Xi7g9WkJVb5HU4xkScZ2ZlUFxWwe0NYGCjYASTLJ4Ukb+IyLEiMrTyEWglVS0DbgSm4/zQv6uqq0XkARGpvK/UI0Az4D0RWeZTvfU+zoWAK4HlwPLKGxuauhEfG8PTlwwlPi6G699YwoFS6x1lTKXNO/fx5oJNnD+sCz3aRvdwqcEKphrqSOAyYDw/V0Op+7pGqjoNmFZl2n0+z0+qZr1y4LdBxGYOQ6dWTZhywWCuenUh93+6mn+cd5TXIRkTEZ74cj0I3HJSb69DiRjBJItzgB5u91cTZcb1a8f1J/Tk2dmZHJ2SxHnDugReyZgolrFjDx8u2czVo1Pp2DL6h0sNVjDVUMuBVqEOxHjntpP7MCI1iXs/XkX69j1eh2OMpx6bkU6T+Fh+d0JPr0OJKMEki/bADyIyvZZdZ009ERcbw1MXDyGxUSzXv7GEvcVlXodkjCdWbN7FF6t+5NdjetCmgQyXGqxgqqH+EvIojOfat2jMkxcNYfJLC7jno5VMuXBwg75a1TRMj0xfR+um8fx6TKrXoUScYK7gtrErGojRvZL5/Yl9mPJlOsf0aMPFI7p5HZIxYTMvM49v1+dx78QjaN6AhksNVk3Dqu7B/xXXAqiqtghZVMYzN47vxaKNBfzl09Uc1aUlAzrZxUgm+qkqD/9vHR1bNmbyyO6BV2iAahopr7mqtvDzaG6JInrFxghTLhxM66bx3PDGEvYcKPU6JGNCbuaa7SzL2cUtJ/ZucMOlBqvhDfdkAkpu1oinLh5Kzs793P3BSlTtTiomepVXKI/OWEeP5ER+ZV3Hq2XJwvg1IjWJP5zSl89XbuP1+Ru9DseYkPl0+RbStxdx2yl9iGuAw6UGy94ZU63fHt+D8f3a8eDna1ies8vrcIypcyVlznCp/Tu24PSBHb0OJ6JZsjDViokRHjt/EO2aN+b6N5ZQuM/aL0x0eWfhJnIK9nPHhL4NdrjUYFmyMDVqnZjA05cMYceeA9z+3jJrvzBRY19JGVO/zmBEShIn9GnrdTgRz5KFCWhIt9b88bQj+HLtDl78NsvrcIypE6/O20DunmLunNDXLkANgiULE5SrRqcwYUAH/vm/dSzaUOB1OMYclsJ9pfxrdibj+7VjeErDHi41WJYsTFBEhIfPP4rOrZpw45tLyS8qDrySMRHq+TmZ7D5Qxh9OseFSg2XJwgStReN4nr10KAX7Srj13eVUVFj7hal/duw5wCtpGzhrUCf6d7Lri4NlycLUysDOLfnLmf2Zk57Ls7MzvA7HmFp75usMSsoruNWGS60VSxam1i4Z0Y2zBnXi8ZnpzMvM8zocE2H2l5Rz6zvL+Nc3mewviazhenMK9vHm95u4YHhXUpMTvQ6nXrFkYWpNRPjbuUeSkpzIzW8tY8eeA16HZCKEqnLvx6v4aOkW/vHFDxz/yCxem7eB4rLISBpTvkwnRoRbTrThUmvLkoU5JM0axfHspUMpKi7llreWUW7tFwZ4d1EOHyzZzM3je/HedcfSIzmRv3y6mvGPfsO7C3MoK6/wLLb07Xv4aOkWrhiVQoeWjT2Lo76yZGEOWb8OLfi/SQOZn5XPk1+mex2O8diarbu575PVjO7VhltO6sPRKUm8fe1I/n3NCJKbJXDnBys4ZcocPl2+1ZPOEY9OX0ezhDh+N9aGSz0UlizMYTl/eFfOH9aFp2Zl8E16rtfhGI/sPlDK9W8spmWTeJ64cAix7q0zRIQxvdvy8Q2jeeGyYcTHxnDzW0s5feq3zFyzPWx3BFi6aScz1mznN8f3oHViQlj2GW0sWZjD9sCkgfRp15xb31nGtsL9XodjwkxVuev9FeTs3M/TlwylbfODx64WEU4Z0IEvbhnDkxcN5kBpOb95fRHnPDuPuevzQp40Hpm+jjaJCVx9nA2XeqgsWZjD1iQhlmcnD6W4tJyb3lxKqYf10ib8XknbwBerfuSOU/syIrXmq6FjYoRJgzsz87ax/OPcI9mx+wCTX1rAxS9+x+KNobkzwNz1eczLzOf6cb1o1ijgSNKmGpYsTJ3o2bYZfzv3SBZt3Mmj09d5HY4JkyWbdvK3aWs56Yh2XDumR9DrxcfGcNGIbnz9hxP4y5n9ydhRxHnPzeeqV75n1ZbCOotPVXlk+g90atmYS4+xMeUPhyULU2cmDe7M5JHdeH5OFjPXbPc6HBNiO/eWcOMbS+jQsjGPnT/4kG7x3Tg+lqtGpzLnznHcOaEvSzbt4oyn5nL9G4vJ2LHnsGOcvno7yzcX8vuT+thwqYfJkoWpU/dO7M+ATi24/d1l5BTs8zocEyIVFcqt7y4jr6iEZy8dSsum8Ye1vaYJcVx/Qi/m3DmOm8f34pt1uZwyZQ63v7v8kL9H5RXKYzPW0bNtIucO7XxY8RlLFqaONY6P5dlLh6IKN765hJIya7+IRs/OzmD2ulz+fGZ/jurSqs6227JJPLed0pc5d47jmuNS+WzFVsY/Npt7P17J9t21u/jzo6VbWL+jiNtP6WvDpdYBewdNneveJpFHzj+K5ZsL+du0tV6HY+rYvMw8Hp+ZzlmDOjE5RO0AbZo14p6J/ZlzxzguGN6Vt7/P4fiHZ/HQ52so2FsScP3isnKmzEznyM4tOW1gh5DE2NBYsjAhMWFgR64encqr8zYwbeU2r8MxdWT77gPc/NZSUpMT+fu5R4Z80KAOLRvz0DlH8vXtJzDxqI68NDebMf/8msdnrGP3geqH+X1rwSa27NrPHafawEZ1xZKFCZm7T+vH4K6tuPP9FWzI2+t1OOYwlZVXcNNbS9lbXM5zk4eRGMZuqN3aNOXxCwYz49bjGdu3LVO/zmDMP2fx7OwM9pWU/WLZvcVlPD0rg5E9khjTOzlsMUY7SxYmZBLiYnjm0qHExQrXv7GEA6WRcTM5c2genZHO99kFPHTOQPq0b+5JDL3aNefZS4fx35uOY2i3Vjz8v3Uc//BsXk3L/ulmha+kZZNXVMIdp/azUkUdsmRhQqpzqyY8fsEg1mzbzV8/W+N1OOYQfbV2O//6JpOLR3Tl3KFdvA6HgZ1b8spVI3j/umPp2TaR+z9bw7hHZvP6/A08PyeLk45ox7Durb0OM6qENFmIyAQRWSciGSJyt5/5t4nIGhFZISJfiUh3n3ndRGSGiKx1l0kJZawmdMb3a8/vTujJW99v4qOlm70Ox9RSTsE+bnt3Of07tuAvZw7wOpxfGO5zs8K2LRpz3yerKSou4w+n2nCpdS1klY4iEgs8A5wMbAYWisinqup7erkUGK6q+0Tkd8DDwIXuvNeBh1R1pog0A6wPZj12+8l9WLxhJ3/6cBUDO7Wkt0fVGKZ2isvKueHNJVRUKM9NHhqRF7ZV3qzwuF7JfLV2B3tLyujXwYZLrWuhLFmMADJUNUtVS4C3gUm+C6jqLFWtvOLmO6ALgIj0B+JUdaa7XJHPcqYeiouN4alLhtA0IZYLnp/P+4s3h+2Oo+bQPfT5WlZsLuSR8wfRvU1kjywnIpzUvz2TBtsFeKEQymTRGcjxeb3ZnVada4Av3Od9gF0i8qGILBWRR9ySyi+IyLUiskhEFuXm2u2xI137Fo15+9qRpCYn8of3lnPJiwvIyi3yOixTjc+Wb+X1+Rv59XGpTLBrFRq8UCYLf90Q/J5KishkYDjwiDspDhgD/AE4GugBXHnQxlRfUNXhqjq8bdu2dRGzCbHe7Zvz/nWjePDsgazaWsiEJ79l6lfrI2bYTePIzC3i7g9WMKx7a+46rZ/X4ZgIEMpksRno6vO6C7C16kIichJwD3CWqhb7rLvUrcIqAz4GhoYwVhNGMTHC5JHd+eq2sZzcvz2Pz0xn4tS5fJ8dmltUm9rZX1LO9f9ZQqP4WJ6+ZAjxdqsMQ2iTxUKgt4ikikgCcBHwqe8CIjIEeB4nUeyosm5rEaksLowHrN9llGnXojHPXDKUV648mv0l5Vzw/Hzu/mAFu/YFvp2DCQ1V5d6PV5G+Yw9PXDiYji2beB2SiRAhSxZuieBGYDqwFnhXVVeLyAMicpa72CNAM+A9EVkmIp+665bjVEF9JSIrcaq0XgxVrMZb4/q1Y+Ztx/Pb43vw3uLNnPjYN3y8dIs1gHvg3UU5fLBkMzeN783xfaxq1/xMouUfcvjw4bpo0SKvwzCHafXWQv700SqW5+xiTO9kHjx7YMT3wokWq7cWcs6z8xiRksRrV4/4aRxtE91EZLGqDg+0nFVGmogyoFNLPvzdKP561gCWbtrFKVPm8MysDBuqNcR2Hyjl+jeW0LppPE9cNNgShTmIJQsTcWJjhCtGpfDlbWMZ17cdj0xfxxlT54ZsjOZQ27WvJKLH9VBV7nxvBZt37ufpS4aS3KyR1yGZCGSjl5uI1aFlY/512TBmrtnOfZ+s4rzn5nP24E6M69eOkT3a0L5FY69DPEhFhZKZW8SijTtZtGEnizcWsCF/H8nNGjF5ZDcuPaY7bZtH1o/xy2kb+N/qH/nT6f04OiXJ63BMhLI2C1MvFBWXMWVmOu8szKGo2LkldWpyIsekJjGyRxuO6ZHkSc+dA6XlrNhcyMINBSzeuJMlm3aya58zzkKbxASGdW/NoK6tWLShgFnrckmIjWHS4E5cNTqV/p28vyXF4o07ufD5+Yzr144XLhtmd2ltgIJts7BkYeqVsvIK1mzbzYKsAr7Lyuf7DQXsOeAkj+5tmvokjzZ0blX3ySOvqPinEsOijTtZtaWQ0nLnf6hn20SGd09iWEprhndvTWpy4i9+fDNzi3g1bQPvL97M/tJyRvZI4urRqZx4RHtP2ggK9pYwceq3xMUK/71pDC2bHN442qZ+smRhGoTyCmXttt18l5XPguwCvs8uoHC/c2bfpXUTJ3G4CaRrUlPAqaMvLqtgb3EZe4vL2VtSxt7iMoqKy9hXUk5RsfPa9/mufaWs3FJItjuIU0JsDEd1aekmhiSGdW9NUmJCUDEX7ivlnUWbeG3eRrbs2k+3pKZcOSqF84d3oXnj8PxgV1QoV766kO8y8/ngd6M4skvLsOzXRB5LFqZBqqhQfvhxDwuy81mQVcCC7Hx2utVCSYkJlJVXsK+knLKK4L73cTFCYqM4mjWK44iOLTg6pTXDU1ozsHNLGsUd3h1Yy8ormLFmOy/PzWbRxp00axTHBcO7cuWoFLq1aXpY2w7kqa/W89jMdB48eyCTR3YPvIKJWpYsjMFJHut3FLEgO581W3fTKC6GxEZxPyWApgmxNHNfJzaKdf4m/Pz6cBNCsJbn7OKVtGz+u2Ib5aqcfER7rj4ulWNSk+q8HWFeRh6TX1rAmYM68cSFg62dooGzZGFMPbR99wH+PX8jbyzYyM59pfTv2IKrj0vlzEEd6yRxbd99gIlTv6VV0wQ+uWF0WMfRNpHJkoUx9diB0nI+XrqFl9OySd9eRHKzBCaP7H5YXW/Lyiu45MUFrNxSyKc3jrYBqAxgycKYqKCqpGXk83JaNl//sIOE2BjOGtyJq0anMKBT7Rql//7FWp7/JosnLhzM2UNsgCDjCDZZWBnUmAgmIhzXO5njeieTlVvEq/M28N6izby/eHOtut7OXLOd57/J4pJjulmiMIfEShbG1DO17XqbU7CPiVO/pVubprx/3aiIHEfbeMeqoYyJcsF0vS0uK+f8f80nO28vn980JuRdck39Y9VQxkS5uNgYTj+yI6cf2fGnrrevz9/AK/OyOemI9lw9OpVpK7exYnMhz182zBKFOSxWsjAmilTtegvwmzGp3DOxv8eRmUhl1VDGNGCVXW+z8/byh1P72jjaplpWDWVMA9Y4PpaLRnTzOgwTRex0wxhjTECWLIwxxgRkycIYY0xAliyMMcYEZMnCGGNMQJYsjDHGBGTJwhhjTECWLIwxxgQUNVdwi0gusNHrOA5BMpDndRBhZsfcMNgx1w/dVbVtoIWiJlnUVyKyKJhL7aOJHXPDYMccXawayhhjTECWLIwxxgRkycJ7L3gdgAfsmBsGO+YoYm0WxhhjArKShTHGmIAsWRhjjAnIkkUYicitIrJaRFaJyFsi0lhEUkVkgYisF5F3RCTB6zjrkojc4h7vahH5vTstSURmusc8U0Raex3n4RKRl0Vkh4is8pnm9zjFMVVEMkRkhYgM9S7yQ1PN8Z7vfs4VIjK8yvJ/dI93nYicGv6ID181x/yIiPzgfo4fiUgrn3n1/ph9WbIIExHpDNwMDFfVgUAscBHwT2CKqvYGdgLXeBdl3RKRgcBvgBHAIOAMEekN3A185R7zV+7r+u5VYEKVadUd52lAb/dxLfBcmGKsS69y8PGuAs4F5vhOFJH+ON/1Ae46z4pIbBhirGuvcvAxzwQGqupRQDrwR4iqY/6JJYvwigOaiEgc0BTYBowH3nfnvwac7VFsoXAE8J2q7lPVMuAb4BxgEs6xQpQcs6rOAQqqTK7uOCcBr6vjO6CViHQMT6R1w9/xqupaVV3nZ/FJwNuqWqyq2UAGzglEvVLNMc9wv9sA3wFd3OdRccy+LFmEiapuAR4FNuEkiUJgMbDL58u2GejsTYQhsQo4XkTaiEhT4HSgK9BeVbcBuH/beRhjKFV3nJ2BHJ/lou1zr6qhHO/VwBfu86g7ZksWYeLWV08CUoFOQCJOdURVUdOXWVXX4lSzzQT+BywHympcqWEQP9Oi5nP3I+qPV0Tuwfluv1E5yc9i9fqYLVmEz0lAtqrmqmop8CEwCqcKIs5dpguw1asAQ0FVX1LVoap6PE4Rfj2wvbLaxf27w8sYQ6i649yMU8KqFHWfexVRfbwicgVwBnCp/nzhWtQdsyWL8NkEjBSRpiIiwInAGmAW8Ct3mSuATzyKLyREpJ37txtO4+dbwKc4xwpReMw+qjvOT4HL3V5RI4HCyuqqKPUpcJGINBKRVJyG/e89jqlOiMgE4C7gLFXd5zMr+o5ZVe0RpgfwV+AHnLr8fwONgB44X6IM4D2gkddx1vExf4uTFJcDJ7rT2uD0Dlrv/k3yOs46OM63cNqiSnHOKq+p7jhxqiieATKBlTg95Dw/hjo43nPc58XAdmC6z/L3uMe7DjjN6/jr8JgzcNomlrmPf0XTMfs+7HYfxhhjArJqKGOMMQFZsjDGGBOQJQtjjDEBWbIwxhgTkCULY4wxAVmyMA2KiNzj3hl1hYgsE5Fj3Om/d29JUt16/8+9ORwiUnQY+98gIitFZLmIzBCRDoH2Z0wksK6zpsEQkWOBx4ETVLVYRJKBBFXdKiIbcK53yPOzXqyqlvu8LlLVZocYw0/7EZG/Ac1U9eaa9mdMJLCShWlIOgJ5qloMoKp5bqK4Ged+XbNEZBY4CUFEHhCRBcCxIjLbzxgNySIyX0Qmuq/vEJGFbqnlr0HEMwfoFWh/IjJBRJa4pZGv3GmJ7vgKC0VkqYhMcqcPEJHv3VLTCveW8MYcNksWpiGZAXQVkXQReVZExgKo6lSc+/aMU9Vx7rKJwCpVPUZV51bdkIi0Bz4H7lPVz0XkFJxbOowABgPDROT4APGcgXMFd7X7E5G2wIvAeao6CDjfnXUP8LWqHg2MAx4RkUTgOuBJVR0MDMe50tiYw2bJwjQYqloEDMMZcCgXeEdErqxm8XLgg2rmxePcvuNOVZ3pTjvFfSwFlgD9cJKHP7NEZBnQAvh7gP2NBOaoMyYCqlo5nsIpwN3udmYDjYFuwHzgTyJyF9BdVfdXE4MxtRIXeBFjoofbFjAbmC0iK3Fu8Peqn0UP1NBuUIYzFsmpOAM6gXO/p7+r6vNBhDHOT9tIdfsT/N/aWnBKG1UHG1rrVmVNBKaLyK9V9esgYjKmRlayMA2GiPStUoc/GNjoPt8DNA9yU4oz0E0/EakcKnU6cLWINHP31bnyjruHaT4w1r1zKSKS5LO/m9w7GCMiQ9y/PYAst2rtU+CoOojBGCtZmAalGfCUiLTCKR1k4FRJAbwAfCEi23zaLaqlquUichHwmYjsVtVnReQIYL77+10ETOYwx+pQ1VwRuRb4UERi3O2dDPwf8ASwwk0YG3DaQC4EJotIKfAj8MDh7N+YStZ11hhjTEBWDWWMMSYgSxbGGGMCsmRhjDEmIEsWxhhjArJkYYwxJiBLFsYYYwKyZGGMMSag/w/l0we2rVQuBAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(strikes, impvols_put)\n",
    "plt.xlabel('Strike Prices')\n",
    "plt.ylabel('Implied Volatilities')\n",
    "plt.title('AAPL Put Implied Volatilities expiring in 78 days')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
